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Executive  Summary 


Recent  conflicts  and  peacekeeping  missions  have  revealed  the  need  for  novel  warheads 
able  to  selectively  release  the  blast  and  fragmentation  energy  against  prescribed  targets  to 
minimise  collateral  damage.  Analysis  of  these  new  warheads  is  impossible  without 
consideration  of  multi-phase  flows  involved  in  the  warhead  detonation.  Warhead 
components,  can  include  composite  and  inhomogeneous  explosives,  and  fragmentation 
warheads  include  fragment  particles,  which  could  be  treated  as  an  additional  phase  of  the 
warhead  flow.  It  also  provides  the  potential  capability  of  tailoring  the  blast  and  impact 
energy  release  at  a  given  time  and  location  that  would  be  of  significant  benefit  to  the 
optimisation  of  the  warhead  effects  against  specified  targets. 

Numerical  analysis  of  the  warhead  effects  involves  at  least  two  stages: 

i)  construction  of  kinetic  relationships,  which  are  responsible  for  the  internal  and 
chemical  processes  in  the  products  involved  in  the  multi-phase  flow.  The 
fitting  of  these  relationships  to  available  test  data  also  involves  their 
verification  with  calculations,  employing  an  accurate  numerical  scheme;  and 

ii)  a  hydrocode/CFD  study  of  the  target  effects  with  a  multi-dimensional  code 
(likely  to  be  a  commercial  code  with  the  possibility  of  incorporation  of  the 
model),  which  employs  the  multi-phase  model  with  the  closing  kinetic 
relations  verified  during  the  previous  stage. 

The  verification  process  is  likely  to  be  conducted  with  a  one-dimensional  code  because 
the  space  accuracy  needed  for  the  verification  is  hard  to  achieve  in  the  multi¬ 
dimensional  case  due  to  resource  limitations.  For  analysis  of  the  processes  involved  in 
the  warhead  detonation  we  need  a  numerical  method,  which  would  describe 
accurately  elementary  processes  (shock  and  rarefaction)  composing  the  whole  picture 
of  expansion  of  the  detonation  products  and  their  interaction  with  &e  target.  The  only 
scheme,  which  explicitly  involves  solutions  adjusted  to  the  elementary  problems 
(Riemann  problems)  is  the  Godunov  method  along  with  the  family  of  follow-up 
schemes  (TVD,  B.  van-Leer  scheme,  etc).  To  apply  this  scheme  to  the  equations  of  the 
model,  certain  requirements  are  necessary  for  the  model's  system  to  be  satisfied;  the 
equations  should  be  written  out  in  the  form  of  conservation  laws  and  elementary 
solutions  should  be  designed  when  building  up  an  appropriate  Riemann  solver. 

Unfortimately,  among  many  models,  having  been  developed  recently,  only  a  few  can  be 
formulated  in  the  form  of  conservation  laws.  Very  significant  progress  in  this  direction  has 
been  made  by  E.  Romensky,  who  proposed  the  conservation  law  formulations  for  a  large 


variety  of  models,  including  a  two-phase  model  from  which  the  present  consideration 
starts. 

The  present  report  considers  the  modelling  of  two-phase  flows  and  suggests  a  new 
formulation  of  the  model  resulting  in  a  thermodynamic  identity  that  is  applicable  to  the 
case  of  one-dimensional  flows.  This  formulation  establishes  a  clear  link  between  the 
pressure  and  energy  definitions,  embracing  the  diffusive  constituents,  which  are  widely 
used  in  the  theory  of  mixtures  [3],  through  this  thermodynamic  identity.  The  present 
formulation  is  more  convenient  for  construction  of  a  Riemann  solver  and  its  use  in  the 
Godunov  scheme.  This  formulation  is  extendable  to  the  case  of  multiple  phases  with 
complete  thermodynamic  closure  of  the  model,  using  only  thermodynamic  potential. 

Results  of  the  present  study  are  important  for  construction  of  the  algorithms  and 
numerical  methods,  which  are  necessary  for  the  verification  of  kinetic  equations  involved 
in  the  process  of  detonation  of  advanced  warheads.  Extension  of  the  model  to  the  case  of 
multiple  phases  is  critical  for  analysis  of  real-life  warheads,  because  novel  volumetric  and 
fragmentation  warheads  involve,  as  a  rule,  more  than  two  phases  in  actual  engagement 
scenarios. 
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1.  Introduction 


Many  modem  warheads  employ  single-  or  multi-stage  initiation  designs  involving  several 
energetic  and  filling  constituents.  Thus,  chemical  and  mechanical  inhomogeneities  are 
widely  used  for  arranging  a  tailored  energy  release  such  as  delayed  reaction/  initiation, 
afterbum,  etc.  Such  the  complex  energetic  materials  involve  multiple  reactive  components 
and  require  description  with  models,  which  are  capable  of  calculating  multi-phase  flows. 
Interest  in  multi-phase  flows  is  spread  world-wide  due  to  application  of  CFD  modelling  to 
processes  of  combustion,  heterogeneous  detonation,  and  to  processes  in  gas-liquid  and 
bubble-liquid  mixtures.  To  simulate  such  the  processes,  models,  comprising  of  the 
conservation  laws  for  mass,  momentum  and  energy,  associated  with  each  of  the  phases, 
are  extremely  popular;  the  conservation  laws  are  formulated  for  partial  characteristics 
(additive  characteristics  with  respect  to  the  common  volume,  containing  several  phases) 
and  interconnected  by  the  exchange  terms.  However,  this  approach  is  not  very  convenient 
because  it  involves  description  of  every  phase  that  multiplies  the  number  of  equations  by 
the  factor  equal  to  the  number  of  phases.  On  the  other  side,  incorporation  of  such  models 
in  a  hydrocode  is  hard  because  the  majority  of  commercial  hydrocodes  operate  with  a 
single  system  of  conservation  laws  complemented  with  so-called  constitutive  equations  (in 
the  CFD/ engineering  chemistry  communities  they  are  usually  called  kinetic  equations). 
Many  attempts  to  consider  multi-phase  medium  as  an  averaged  one  have  been  made, 
including  a  classic  monograph  by  Truesdell  [3].  However,  a  closed  thermodynamic 
formulation,  resulting  in  an  efficient  practical  realisation,  has  not  been  proposed  at  that 
time.  A  variety  of  models  have  been  recently  developed  in  several  papers  [4, 5].  However, 
they  are  not  in  the  form  of  conservation  laws,  that  complicates  analysis  of  their 
thermodynamical  correctness  and  makes  application  of  the  Godunov  scheme  difficult. 

The  present  work  employs  a  consistent  approach,  enabling  us  to  derive  equations  in  the 
form  of  conservation  laws;  one  of  the  first  realizations  of  this  approach  has  been  published 
in  [6].  It  invokes  the  mass  averaging  over  two  phases,  so  the  effective  averaging  parameter 
involved  is  the  mass  fraction  of  one  of  the  two  phases.  Realization  of  this  approach  as  a 
computer  code  [6]  resulted  in  significant  numerical  difficulties  associated  with  the  phase 
exchange  (convection)  in  the  areas  of  high  gradients  between  the  phases.  A  note  by 
Drumheller  [7]  was  important  for  understanding  that  one  more  parameter  associated  with 
the  phase  concentration  should  be  involved  in  the  processes,  which  accompany  phase 
exchange,  reaction  processes  between  the  phases,  etc.  This  resulted  in  introduction  of  both 
mass  and  volumetric  concentrations  in  the  averaging  process  and,  as  a  result,  several 
successful  models  [2,  8]  based  on  this  approach  have  appeared:  the  model  [2]  is  an 
extension  of  the  model  [6]  for  the  case  of  two-phase  media  with  the  velocity 
nonequilibrium  (drag)  between  the  phases  and  [8]  is  a  single-velocity  two-phase  model 
with  the  temperature  nonequilibrium  resulting  in  a  complete  thermodynamic  identity, 
enabling  one  to  formulate  the  Gibbs  potential  in  its  classical  form.  Formulation  [2] 
provides  a  thermodynamically  correct  model,  but  known  solutions  of  the  Riemann 
problem  cannot  be  easy  applied  because  the  pressure  and  energy  characteristics  involved 
in  the  jump  conditions  are  not  deduced  from  a  single  potential.  The  generalized  pressure 
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and  energy  can  be  introduced  [3],  which  involve  the  diffusion  fluxes;  however,  these  force- 
energy  characteristics  have  not  been  linked  within  the  formulation  [2]. 

The  present  publication  is  an  enhanced  formulation  of  the  model  [2]  for  the  one¬ 
dimensional  case,  allowing  us  to  formulate  a  complete  thermodynamic  identity,  jump 
conditions  in  a  convenient  form  for  application  to  the  Riemann  problem,  and  to  generalize 
the  model  for  flue  case  of  multiple  phases,  linking  directly  the  presentation  [3]  for  energy 
and  pressure  with  involvement  of  diffusive  components.  The  multi-phase  generalization  is 
particularly  important  for  many  applications  because  a  typical  warhead,  for  instance,  of 
volumetric  action  may  involve  at  least  three  phases:  gaseous,  dispersed,  and  liquid 
products. 


2.  Model  of  two-phase  flows 


Let  us  denote  the  average  density  of  a  two-phase  medium  by  p=mfV,  here  m  is  mass  of  a 
representative  volume  and  V  is  the  volume  quantity.  Similarly,  we  can  define  specific 
densities  of  the  phases pi=mil  Vi  and  p2=mil  Va.  Multi-phase  theories  usually  deal  with  so- 
called  partial  densities,  which  relate  the  phase  masses  to  the  whole  volume  such  as: 
p\=mi /  V  and  p' 2=1112/  V.  The  partial  characteristics  are  important  because  the  conservation 
laws  for  each  phase  can  actually  be  formulated  only  for  these  characteristics.  For  the  case 
of  media  with  phases,  which  are  capable  of  an  exchange  with  mass  and  momentum,  the 
conservation  laws  in  one-dimensional  case  take  the  following  form  for  the  first  phase: 


dp[  dp[u, 

■  “r  ' 

dt 

dp[u^ 


dx 


m 


0  ’ 


+ 


d(p[u^  +p[) 


-  n. 


0  ’ 


dt  dx 

+  u 


dt 


■  +  ■ 


dx 


=  0, 


(1) 


and  for  the  second  phase: 


dp2  ^  dp'jUi 
dt  dx 


-Wn 


dp'^u^  .  sjpWi  +p'2)_ 


dt  dx 

+  u  2/2)  ,  ^[^>2(^2  +  u  ll2)  +  p'2U^]_ 
dt  dx 


=  0. 


(2) 


Here  m*o  is  the  mass  exchange  rate,  n‘o  is  the  momentum  exchange  rate,  Ui  (i=l,2)  are 
velocities  of  the  phases,  p'l  and  p'2  are  partial  pressures  within  the  phases,  ei  and  62  are 
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Specific  internal  energies.  Let  us  denote  T  -  temperature  and  S  -  specific  entropy,  then  the 
thermodynamic  identity 

T  dS  =  de  +  p  d  V  =  de  -p  dp/ p  2 ,  (3) 

being  applied  to  each  of  the  phases,  enables  us  to  calculate  partial  pressure  and 
temperature: 


if  a  dependence  of  specific  energy  on  p'  and  S  is  given: 

=  ei(p'i,Si)  ,  62  =  eiip' ifSi)  .  (5) 

It  should  be  noted  that  definition  of  the  partial  pressure  is  based  on  the  application  of  the 
thermodynamic  identity  with  respect  to  the  partial  density.  Thus,  the  traditional  approach 
to  calculation  of  two-phase  flows  is  to  calculate  the  systems  (1)  and  (2),  pre-selecting  the 
exchange  terms  m'o  and  n’o,  and  tabulating  the  'equations  of  state'  in  the  form  (5)  (for  the 
sake  of  convenience,  we  call  the  relations  like  (5)  as  equations  of  state),  using  (4)  for 
calculation  of  pressure  and  temperature. 

The  procedure  of  averaging,  having  been  employed  in  [6],  involves  introduction  of 
averaged  density,  pressure,  and  velocity.  On  the  present  stage  we  do  not  individualize  the 
thermal  characteristics  of  the  phases;  an  example  how  it  could  be  done  for  a  single-velocity 
material  was  shown  in  [8];  therefore,  we  consider  specific  entropy  to  be  common  for  the 
both  phases.  We  introduce  [1,  2,  8]  the  mass  concentration  of  the  first  phase  as  c  =  Ci  = 
nii/ni,  then  for  the  second  phase  C2  =  1112/  m  =  1  -  c.  We  can  also  determine  the  partial 
densities,  because 

p'i=mi/V={nii/m)im/V)=  pc  ,  p'2=p(l-c)  .  (6) 

Specific  energy  is  an  extensive  variable,  therefore,  for  a  volume  containing  both  phases  the 
average  specific  energy  is 

e  =  cei  +  {1 -  c)  62. 

Introducing  volumetric  concentration  of  the  first  phase  as  0i = 0  =  Vi/  V,  we  can  recalculate 
the  specific  densities  of  the  phases 

pi  =  mi/  Vi=  {mi/m)  -{m/  V)  iV /  Vi)  =  pc/6  ,  P2=  p{l  -  c)/(l  -  6)  •  (8) 


3 


DSTO-TR-1510 


Using  (7)  and  (8),  we  can  build  up  the  equation  of  state  for  the  averaged  state,  employing 
the  'local'  equations  of  state  e\  =  ei{pi,  S)  and  €2=  e2(p2,  S): 

e{p,  c,  d,S)  =  c  eiipc/d,  S)  +  (1  -  c)  eiipO-  -  ^)/0-  "  ^)'  ■  (^) 

Averaged  velocity  u  is  introduced  as 

U  =  Clll  +  {1  ~C)  tl2  ,  (10) 

and  the  velocity  difference  between  the  phases,  which  is  proportional  to  the  so-called 
diffusion  velocity  [3],  was  introduced  [2,  6]  in  the  form 


IV  =  Ml  -  M2  .  (11) 

Relations  (10)  and  (11)  allow  us  to  calculate  local  velocities  via  the  averaged  velocity  and 
the  velocity  difference: 

Ml  =  M  +  (1  -  c)  w  ,  U2-11-CW  .  (12) 

From  (6)  and  (10);  p\  +  p'2  =  p  and  p’l  mi  +  p'2 112  =  p  m;  using  these  relations  and  summing 
up  the  continuity  equations  in  (l)-(2),  we  can  obtain  the  continuity  equation  for  the 
averaged  variables: 

^  +  ^  =  0.  (13) 

dt  dx 


Rewriting  the  continuity  equation  of  (1)  with  the  use  of  (6)  and  (11),  we  can  obtain  the 
kinetic  equation  for  the  mass  concentration: 


dpc  d\f>uc  +  pc(l  -  c)w] 
dt  dx 


* 


(14) 


Kinetic  equation  for  the  volume  concentration  is  chosen  in  usual  form  of  conservation 
within  the  liquid  volume  [1,  2,  6]: 


dp6  ^  dpu6  _  ^ 


dt 


dx 


(15) 


here  <l>  is  a  function  responsible  for  the  process  of  phase  compaction. 

For  further  derivations  we  have  to  calculate  pressure  in  the  averaged  medium.  Firstly,  we 
link  the  local  pressures  and  densities  with  the  partial  ones.  In  addition  to  (6),  another 
calculation  of  the  partial  density  for  the  first  phase  gives 
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p\=mi/V  =(mi/Vi)-{Vi/V)=  pid  ,  p'2=p2(l-0)  .  (16) 

We  consider  an  alternative  to  (5)  presentation  of  the  equations  of  state  in  the  form  ei  =  ei(pi, 
S),  ez  =  e^(p2,  S).  Then,  from  (4)  and  (16)  it  follows 


Pi  . 


where  the  following  denotations  for  the  local  pressures  are  used: 

_  /  _  2  de./ 

Pi  -  Pi  ,  Pi -  Pi  . 

with  the  equations  of  state  given  in  the  form: 


(17) 


(18) 


ei-ei(pi,S)  ,  e2-e2(p2,S).  (19) 

The  momentum  conservation  laws  (second  equations  of  ihe  systems  (1)  and  (2))  give  the 
following  momentum  equation  for  the  two-phase  medium: 

dpu  ^  d[pu^+p  +  pc(\-c)\v^]  ^ 

dt  dx 

Here  the  following  relations  have  been  used 


P'l  Ul+  p'2U2  =  pU  , 

P’l  (Mi)2  +  P'2  (1/2)2  =  plP  +  p  C(1  -  c)z//2, 

and 


(21) 


P=p‘l  +  P'2. 

The  latter  can  also  be  obtained  by  direct  differentiating  of  (7)  over  p  and  gives  the  relation 
p=6pi  +  {1-6)p2. 

The  momentum  conservation  laws  of  (1-2)  can  also  be  rewritten  in  the  following  form 
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du.  du,  1  dp[  (  *  » 

^  +  w,^  +  — =  K -m,u,}/p,  , 
dt  dx  Pi  ox 


du. 


■  +  u. 


du. 


■  + 


■V^  =  -("o  -^1^2)1  P'2  ■ 


dt  ^  dx  /?2  Sx 
When  using  the  relation 
(Mi)2  -  (M2)2  =  2uw  +  (1  -  2jc)i{P  , 

we  can  derive  the  equation  for  the  velocity  difference: 


dw  ^ 
dt 


mv 


+  (l-2c)'^y2J^  1  dp'i  1  dp'2 
dx  p[  dx  p2  dx 


(22) 


where 

pc  p(l-c)  ■ 

To  simplify  the  equation,  a  chemical  potential  n  is  introduced  [2],  which  could  be  equal  to 
the  classical  Gibbs  potential  if  the  energy  exchange  is  properly  introduced  [8].  The 
potential  is  defined  as  follows 


U  =  Cc  =  fix  -  C2  +  pi/  pi  -  P2/  P2  • 

Noting  that  pil  pi  =  p' it  p'i,'w^  can  simplify  the  last  kinetic  equation  as  follows: 


div  ^ 
dt 


\uw  + 


dx 


(23) 


The  last  set  of  conservation  laws  in  (l)-(2)  deals  with  the  energy  conservation.  To  derive  it 
for  the  averaged  medium,  we  first  obtain  an  auxiliary  relation: 


c(lll)2  +  (1  -  C)(U2)2  =  Ip  +  c(l  -  C)zy2  . 

Using  this  relation  and  (22),  we  can  calculate  the  following 
p'iui[ei  +  (iXi)2/2]  +  p'lUilei  +  (m2)V2]  = 

=  p{ue  +  tP/2  +  ii  c(l  -  c)io^/2  +  c(l  -  c)io  [ei  -e2+uio  +  0--  2c)iu‘^l2] }  , 


and  from  (21): 
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p'llei  +  {uiy/2]  +  p2[e2  +  (1/2)72]  ^p[e  +  +  c(l  -  c)ivy2]  . 

Finally, 


p'll/l  +  p'2ll2=p\[u  +  (1  -  c)  id]  +  p'llll  -Civ]  = 

=  (p'l+  p'2)  u  +  tw[(l  -  c)  p'l  -cp'2]  =pll  +  lO  [p'i(l  -c)  p'l/  P  \  -p2Cp2l  p'2]  = 

=  pu+piu  [c(l  -  c)  p\/  p'l  -  c(l  -  c)  p'l/  p'2]  =  pu  +  pwc{l-  c)  [pi/  pi  -  p2/  p2]  • 

Let  US  denote 

E  =  e  +  c{l-  c)iv‘^/2,  (24) 

then  the  energy  conservation  equation  takes  the  following  form  for  the  two-phase 
medium: 


8p{e  +  ^^2 )  8\pu{e  +  1+  pw  +  /Oc(l  -  c)ir[u  +  MW  +  (1  -  2cy^ 

dt  dx 


=  0. 


It  is  seen  that  it  is  convenient  to  calculate  the  chemical  potential  from  the  specific  energy  E 
including  the  diffusion  term: 


A  =  Ec  =  ei  -  £2  +  pi/  pi  -  P2/  p2  +  (1  -  2c)  ty2/2  . 

Then  tire  kinetic  equation  for  the  velocity  difference  (23)  takes  the  following  form 

dw  8[uw  +  a]  _  ^ 

8t  8x 


It  can  be  noticed  that  Ew  =  c(l  -  c)iv,  then  the  energy  conservation  takes  the  following  form: 

8p{E  +  ^'y^j^  +  pu  +  pinvE^^  +  pE^A^  ^ 

8t  8x 

Summarizing,  the  complete  system  of  equations  for  a  two-phase  medium  involves 
equations  (13-15),  (20),  and  (25-26): 
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dt  dx 

dpu  ^  d\pu^  +  p  +  pwE,]  _^ 
dt  dx 


ap(£+«^)  +  pu^vE^,  +  pE^^I^ 


dt 

dpc  d[puc  +  pE^,]_  . 

dt  dx 

dp9  ^  dpu6  _  ^ 
dt  dx 

dw  dlmv  +  a] 
dt  dx 


dx 


(27) 


This  is  almost  identically  the  system  of  equations,  describing  a  two-phase  medium,  which 
has  been  obtained  in  [2].  The  system  is  thermodynamically  consistent,  because  if  a 
generalized  equation  of  state  is  given  in  the  form 

E  =  E{p,S,c,  B,  lo) , 

then  the  entropy  evolution  equation  can  be  obtained  from  (27)  and  the  conditions  of 
hyperbolicity  deduced  [2].  However,  within  the  formulation  [2]  there  is  no  a 
thermodynamic  identity  specified,  which  would  clearly  relate  the  forces,  appearing  asp  + 
pwEv,,  to  the  energy  £. 

3.  Thermodynamic  identity 


In  this  section  we  are  reconsidering  the  model  to  derive  thermodynamic  identity  and 
obtain  a  more  convenient  system,  which  could  be  generalized  for  the  case  of  multiple 
phases. 

Let  us  analyze  the  jump  conditions  for  the  system  (27) .  We  neglect  kinetic  rates  (the  right- 
hand  sides  in  (27))  for  derivation  of  the  conditions.  Let  us  denote  the  jump  velocity  by  D, 
and  then  the  jump  conditions  take  the  following  form: 

p  {ii  -  D)  =  const  , 

pii  {ii-D)  +  p  +  pwEw  =  const  , 

p(E  +  tP/2)  {u-D)  +  u{p  +  pivEw)  +  pEwA  =  const  ,  (28) 

pc  {u-D)+  pEw  =  const  , 
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It  is  well  known  that  for  a  two-parametric  medium  (a  medium,  which  can  be  completely 
described  with  two  parameters,  e.g.,  density  and  entropy)  the  conditions  of  continuity  on 
the  contact  jump  appear  as  the  equality  of  pressure  and  velocity.  Here  we  describe  the 
contact  jump  as  a  jump,  moving  with  the  medium,  i.e.,  the  jump  velocity  is  equal  to  tihe 
fluid  velocity.  Let  us  derive  similar  conditions  for  the  present  two-phase  medium.  If  we 
take  ii  =  D,  then  the  contact  conditions  follow  from  (28)  as 

p  +  pxvEw  =  const  ,  pEw  =  const  ,  A  =  const  .  (29) 

It  is  seen  that  the  role  of  pressure  play  the  following  functions:  p  +  pioEw,  pEw,  and  A.  We 
suppose  that  these  functions  relate  directly  to  the  thermodynamic  relations  of  the  medium. 
Let  us  introduce  a  generalized  pressure: 

P-p  +  pioEw  ■  (30) 

It  should  be  noted  that  replacement  of  pressure  by  a  combination  involving  also  the 
diffusive  components  has  been  considered  in  theories  of  mixture  long  ago,  by  Truesdell  [3] 
and  many  other  researchers,  as  well  as  introduction  of  the  diffusion  terms  in  the  internal 
energy  similarly  to  (24).  A  diffusion  force  for  die  second  relation  in  (29)  could  be 
introduced  as 


pEw 


(31) 


Thus,  the  functions,  preserved  through  the  contact  jump,  are  P,  A,  and 

A  quantity,  which  does  not  change  through  any  jump,  is  mass  m  =  p  {ii  -  D).  Then  the 
conditions  (28)  can  be  rewritten  as 

m  =  const  ,  mil  +  P  =  const  ,  m{E  +  ififT)  +  iiP  +  AQ  =  const  , 
mc  +  Q,  =  const  ,  md  =  const  ,  m{w/p)  +  A  =  const  . 

It  is  seen  from  these  relations  that  basic  variables  associated  with  fluxes  of  A,  and  Q  are 
mass  concentration  c  and  iv/p.  Therefore,  it  is  convenient  to  introduce  a  new  variable  (a 
specific  diffusion  velocity): 


(32) 


V  =  w/p  . 

Now  from  (24)  we  can  see  that 

E  =  e(p,  S,  c,  6)+p2c(l-c)uV2  / 


(33) 
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and  the  function  Q  is  calculated  as  =  Ey.  The  function  A  can  be  calculated  as  above  in 
accordance  with  A  =  Ec.  However,  what  is  even  more  important,  the  generalized  pressure 
can  be  calculated  via  E  as  P  =  p^Ep.  The  specific  energy  also  depends  on  6,  so  we  have  to 
introduce  a  function  n  =  Ee.  From  (9)  it  can  be  found 

n  =  c(ei)pi(pc/i9)e  +  (1  -  c)(e2)p2(p(l  -  c)/(l  -  &))q  = 

=  -  (pi)2(ei)pi/p  +  {p2f{e2)pi/p  =  -  (pi  -  p2)/p  . 

Putting  those  equations  together,  we  have: 

T=Es,  P=p2Ep,  A  =  Ec  ,  Q  =  Ev  ,  H  =  Ee  .  (34) 

Thus,  the  following  thermodynamic  identity  takes  place 

TdS  =  dE  +  P  dV  -  A  dc  -  Q  du  -  n  d0  .  (35) 


The  system  (27)  can  be  rewritten  in  the  following  form 


dp  dpu 
dt  dx 


dpu  d  pu^  +  P  _  Q 
dt  dx 


dp\E  + ' 


d\pu\E  +  ^  /^  )+  Pm  +  AQ 


dt 


dx 


=  0 


dpc  d\puc  +  0\  ♦ 

dt  dx 

dp9  ^  dpuO  _  ^ 
dt  dx 


dpv  ^  d[puv^h]  _^^, 
dt  dx 

Concluding,  selection  of  the  potential  in  the  form 


(36) 


E  =  E(p,S,c,  d,  V)  (37) 

along  with  the  identity  (35)  closes  the  model  (36).  In  fact,  the  potential  E  can  be  calculated 
as  earlier,  using  (33)  and  (9)  along  with  the  local  equations  of  state  (19). 
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4.  Model  of  multi-phase  flows 

In  the  present  section  we  shall  follow  more  traditional  denotations,  being  used  in  the 
mixture  theories  (e.g.,  [3]),  which  are  associated  with  the  centre  mass  velocity  and 
diffusion  velocities.  This  choice  is  sometimes  more  convenient  because  it  is  not  associated 
with  a  specific  component  of  a  mixture  and  treats  the  mixture  components  equally.  Thus, 
for  a  n-phase  medium  we  denote  d  as  mass  concentration  of  i-th  component  of  the  mixture 
(i=l, ...  ,n).  We  suppose  that  the  multi-phase  medium  has  n  components  and  summation 
from  i=l  up  to  i=  n  will  be  denoted  by  sign  Z.  As  usual,  d  =  nii/  m,  where  nii  is  the  mass 
fraction  of  the  i-th  component  in  die  volume  V  of  mass  m.  Similarly,  0i  =  Vi/  V  that  gives  Z 
Ci  =  1  and  Z  0i  =  1.  At  the  moment  we  ignore  the  interdependence  of  the  variable  sets  d  and 
0i  at  i=l, ...  ,n;  we  will  recall  this  fact  at  the  end  of  the  section  when  this  interdependence 
is  relevant.  We  introduce  average  velocity  u  exactiy  how  we  have  done  it  in  the  previous 
section: 

u  =  Z  Ci  z/i  .  (38) 

However,  to  preserve  universality,  we  define  the  velocity  nonequilibrium  in  a  different 
manner  in  accordance  with  the  traditional  choice  of  diffusion  velocities,  which  is 
proportional  to  the  well-known  diffusion  fluxes  pi{iii  -  ti): 

X0i  =  Ui-u  .  (39) 

It  can  be  noted  that,  comparing  this  choice  with  the  choice  in  the  previous  section,  the 
diffusion  velocities  would  be  iu\  =  u\  -  u  =  (1-  c)  iv  and  W2  -  ui  -  u  =  -  c  iv.  It  follows 
cwi+{l-c)  102  =  0,  or,  generally: 


Z  Ci  Wi  =  0  .  (40) 

The  diffusion  velocity  interdependence  will  also  be  ignored  at  the  moment.  Keeping  in 
mind  the  results  of  the  preceding  section,  we  select  the  following  variables  as  the  basic 
ones: 


p,  S,  u,  Ci,  di,  Vi.  (41) 

Again,  the  variables  are  independent  at  z=l,...,n-l,  but  at  the  moment  we  consider  the 
overdetermined  set  (41)  at  i=l,...,n.  The  local  densities  are,  similar  to  (8): 


pi  =  pd/di  .  (42) 

From  (39)  iii  =  ii  +  pvi,  that  allows  us  to  determine  the  kinetic  equations  for  Ci: 


^££l+ 

dt 


d{puCj  +  Q, ) 

8x 


* 


(43) 
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here  =  p^CiVi  =  pcm  (Z  =  0  from  (40))  and  m'i  is  the  mass  exchange  rate  such  that  Zm*i  = 
0,  when  assuming  that  for  the  i-th  phase  the  equations  are  chosen  in  the  form  identical  to 
(1)  replacing  the  subscript  index  '1'  by  'i'.  Summing  up  the  continuity  equations  of  (1),  we 
have  the  continuity  equation  for  the  mixture  coincident  with  (13): 

^  +  ^  =  0,  (44) 

8t  8x 

The  volumetric  concentrations  behave  similarly,  regardless  of  the  number  of  phases: 


+  ^  =  .  (45) 

8t  8x 

Let  us  define  a  generalized  internal  energy  of  i-th  phase,  invoking  the  diffusion  energy: 

E\  =  ei  +  w^/  2  =  Ci  +  pk)?-l1  .  (46) 

The  mass  averaging  of  (46)  gives  the  energy  of  the  mixture,  similarly  to  (33): 

E  =  S  CiEi  =  e  +  Z  CiW■?|^  =  e  +  fZ  dv?-!!  ,  (47) 

here  g  =  Z  c^i.  Differentiating  E  over  p,  we  obtain  generalized  pressure  similarly  to  (30): 

P  =  p^Ep  =  p'^p  +  p3Z  CiVp-  =  p'^Cp  +  p  Z  Ciiop-  =  p  +  pZ  drop-  ,  (48) 

where  pressure  can  be  treated  as  in  the  preceding  section:  p  =  Z  p'i  =  Z  6,  pi.  For  derivation 
of  the  momentum  and  energy  equations  for  the  mixture  we  need  several  auxiliary 
relations: 


Z  dup  =u'^  +  'L  dxop , 


Z  pd  (Ci  +  up/2)  =p{e  +  iP  /2  +  Z  diop  /2) , 

Zuip'i  =  up  +  pZdmpi/pi  / 

Z  pdUi  {d  +  up/ 2)  =  pu  (e  +  iP  /2  +  Z  diop  /2)  +  pu  Z  dwp  +  p  Z  d^v^ei  +  p  Z  diop  /2  . 


(49) 


Summing  up  the  momentum  equations  in  the  form  (1),  we  obtain  an  equation  similar  to 
(20)  with  the  generalized  pressure  from  (48): 


8pu  8{pu^  +p)  _  Q 
8t  8x 


(50) 
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When  summing  up,  we  have  to  use  the  momentum  balance  in  the  form  Zn’i  =  0.  To  derive 
the  kinetic  equation  for  the  diffusion  velocity,  we  rewrite  the  momentum  equation  in  the 
following  form 


du.  1  dp] 
+  u, — -+  .  ■ 

dx  /?,  dx 


{n]-m]u)lp] , 


and  the  equation  (50)  in  the  similar  form 


dw  ^  ^  5m  ^  1  dP  _ 
dt  dx  p  dx 


(51) 


(52) 


Using  d  (Ci  +  Pi/ pi)  =  dp'i/p'i,  d  (£  +  P/p)  =  dP/p  and  an  analogue  of  (22)  in  the  form  Ui^  - 
=  2m0i  +  xo?,  we  can  derive  the  required  equation  by  subtracting  (51)  from  (52): 

dxv.  d{mv.+wfj2)  df  /?, 

dt  dx  5x1^  /9,  p^ 


(n]  -m]u,)/p]=W. 


Let  us  introduce  a  chemical  potential  as  follows 


A,. 


£l 

Pi 


(53) 


For  future  references  a  part  of  the  potential,  which  is  not  related  to  a  specific  phase  will  be 
denoted  by  Ao  =  E  +  P/p.  Then  the  equation  for  specific  diffusion  velocity  takes  its  final 
form 


5/7V,  ^  a(pMV,  +A,) 
dt  dx 


(54) 


The  effect  of  change  in  volumetric  concentration  onto  the  specific  energy  is  taken  into 
account  by  introduction  of  Tli,  similarly  to  that  in  the  second  section: 

rii  =  Eei  =  Ci(ei)0i  =  Ci(Ci)pi(pi)9i  Ci{c^pi{p  Ci  / Q'?)  Pi/ p  ■ 


Finally,  the  energy  equation  is  derived  by  summation  of  the  energy  conservation  laws  in 
the  form  (1)  and  use  of  (49): 


dp\^  +  )  5  pm(e:  +  )+  Pm  +  pLc,  w,.  {e.  +  wf  /2  +  p,  /  p, ) 

dt  dx 


=  0. 
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One  more  relation  is  necessary  for  finalizing  the  derivation  of  the  energy  conservation 
equation: 


p  L  cm  {a  +  IV?/ 2  +  Pi/ Pi)  =  Z  pcm  (d  +  iv?/2  +  pi/pi  -  Ao)  =  Z  Ai  , 


here  we  used  definitions  of  ^2iand  Ai,  and  Z  =  0.  Finally,  the  energy  equation  takes  the 
following  final  form 


dp[E  +  “%)  d  pu[e  +  )+  Pu  +  ZA  .Q 

dt  dx 


=  0. 


(56) 


The  complete  system  of  equations,  generalising  (36),  combines  (43-45),  (50),  (54),  (56)  into 
the  following 


dp  dpu 
dt  dx 

dpu  ^  d{pu^+p)  ^^ 
dt  dx 


dp(E  +  d  pu(e  +  Pu  +  ZA,Q, 

dt  dx 


0, 


dpc,  ^  d{puc,+Q,) 
dt  dx 


dp^i  , 
dt  dx 


5 


(57) 


dpv,  ^  d(puv, 
dt  dx 

Concluding  the  section,  we  can  write  out  the  thermodynamic  identity  for  the  multi-phase 
mixture: 


T  dS  =  dE  +  P  dy  -  Z  (Ai  +  Ao)dci  -  Z  t2i  d^i  -  Z  Tli  dOi  ,  (58) 

which  provides  us  with  the  closure  relations: 

T=Es,  P=p2Ep,  Ai  +  Ao=  Ed  ,  ^  =  Evi  ,  Hi  =  Em  .  (59) 

For  completeness  we  have  to  calculate  the  function  Ao  in  the  following  way 
Z  CiEd=  Z  Ci(ei  +  IV?/ 2  +  pi/ pi)  =  Z  Ci(ei  +  iv?/2  +  dipi/{Cip))  =  E  +  p/p  = 

=  E  +  (P  -p  I.  CiW?)/p  =  Ao  -  p^Z  CiV?  , 
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that  is 


Ao=Zci(£ci+pW)  .  (60) 

Regarding  the  entry  of  the  function  Ao  in  the  identity  (58),  we  notice  tiiat  the  summation  is 
conducted  from  1  up  to  n,  so  one  of  the  variables  is  dependent.  For  example,  we  can  select 
the  n-th  mass  concentration  variable  Cn  to  be  expressed  via  the  preceding  ones,  such  as 

Cn  ~  1  “  Cl  -  . . .  —  Cn-1  •  (61) 

It  should  be  noted  that  any  of  these  n  variables  can  be  selected  dependent  on  other  n-1 
ones  simply  by  renumbering  and  assigning  the  number  n  to  the  pre-selected  dependent 
one.  Then  Ae  first  sum  in  (58)  can  be  rewritten  in  the  form,  containing  only  the  first  n-1 
independent  variables: 


Z(Ai  +  Ao)dCi  =  Z'(Ai-An)  ,  (62) 

where  Z'  is  the  summation  sign  from  1  up  to  n-1.  Thus,  the  relation  eventually  allows  us 
to  exclude  the  function  Ao  from  the  consideration. 

Finally,  we  shall  show  that  the  closed  system,  containing  only  n-1  independent  variables, 
is  also  thermodynamically  consistent.  To  do  this,  we  select,  as  for  the  derivation  of  (61), 
that  the  n-th  variables  Cn,  6n,  and  v„  are  dependent  on  the  rest  of  the  set  (41) .  Namely,  from 
(40),  (61),  and  from  the  similar  relation  for  9  we  have 

Cn=l-Z’Ci  ,  en=l-Z'0i  ,  CnVn= -'L'CiVi  .  (63) 

Scalar  functions  in  the  identity  (58)  are  not  affected  by  the  dependencies  (63),  so  the 
relations 


T  =  Es  ,  P=p2Ep 

are  intact.  For  the  next  terms  in  the  identity  we  are  conducting  separate  analyses.  The  first 
one  concerns  Z  (Ai  +  Ao)dci,  which  is  reduced  to  Z'  (Ai-  An ),  according  to  (62).  According 
to  (63),  increment  in  c,  is  also  involved  in  the  change  of  Vn-  Therefore,  for  this  analysis  we 
have  to  expand  the  term  Z  Qi  d^i  as  well: 


Z  Qi  dvi  =  Z'  Qi  dvi  +  £2n  d^n  . 


(64) 


From  (63)  d(Cnyn)  =  -  Z'd  (cii;i)  and  dcn  =  -  Z'dci  that  give 
l^ndCn  +  Cndt^n  =  CndVn  -  l^nZ'dCi  =  -  Z'l'idCi  -  Z' CidVi 


and 
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Cndl^n  =  -  yn)dCi  -  Z'CidUi 


(65) 


From  (63)  and  (64)  it  follows 


Z  Qi  dui  =  Z'  Qi  dfi  -  (Qn/cn)  [2'(yi-  yn)dci  +  Z'cidvi  ]  = 


=  Z’[Qi  -  (Qn/Cn)Ci]  d^i  -  (f2n/Cn)  2'(^i"  yn)dCi 


(66) 


Thus,  according  to  (62)  and  (66),  the  term,  associated  with  the  change  in  Cj  in  the  identity 
(58)  is  actually 


Ai- An- {Qn/Cn)  (Vi-Vn)  .  (67) 

We  have  to  check  out  if  this  term  is  coincident  with  Ed.  Turning  to  (42)  and  (63),  c,  is 
involved  in  pi,  Cn  and  Vn,  differentiating  (47)  over  d  with  this  in  mind  and  using  (65),  we 
have 


£ci  ~  Ei  Ci(ei)pi(pi)d  +  En(CTi)d  "t  Cn(ei)pn(Fn)fT>(Cn)ri  "t  Cn(En)vn(Ei)vi  /Cn 


—  Ei  "t  pi(ei)pi  —  En  “  Pn(^n)pn  —  (^i  ^n)/ Cn  ~  (f^n/Cn)  (^^i  t^n)  , 

which  is  actually  coincident  with  (67).  Next  stage  involves  analysis  of  the  term  with 
increment  in  Vi;  Ais  increment  appears  only  in  the  term  (66)  and  is  equal  to  Qi  -  (Qn/ c„)ci. 
We  have  to  check  out  if  it  is  coincident  with  £«.  The  variable  Vi  is  involved  only  in  Vn; 
differentiating  (47)  over  vi  and  using  (65),  we  can  obtain 


Evi  ~  (Ei)vi  "t  (En)vn(C^n)vi  —  Qn(Ci/Cn) 

that  is  identical  to  -  (Qn/cn)ci.  The  last  analysis  involves  the  term  of  (58)  with  d0i;  taking 
into  account  the  interdependence  of  6i  in  (63),  the  last  term  (58)  is  transformed  as  follows 


Z  Hi  d0i  =  Z '  Hi  d0i  +  Nn  d0n  =  Z'  (Hi  -  nn)d0i  . 


For  the  set  of  independent  variables  0i  we  have  to  find  out  if  Hi  -  fin  is  equal  to  Eei.  To  do 
this,  we  again  differentiate  (47)  over  0i: 

E0i  =  Ci(Ei)gi  +  Cn(En)6n(dn)di  ~  Ci(d)di  ”  Cn(en)9n  “  fli  —  Fin  . 


This  concludes  the  analysis  and  proves  that  for  the  set  of  independent  variables  (41)  at 
i=l,...,n-l  the  identity  can  be  written  out  in  the  following  form 

TdS  =  dE  +  Pdy- 


-  Z '  [  Ai  -  An  -  (Qn/  Cn)  (Vi-  Un)]  dCi  - 


(68) 
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-  Z'  -  Qn  (ci/cn)]  dvi  -  2'  (Hi  -  Hn)  d0i 

Keeping  the  choice  of  the  set  of  independent  variables,  the  next  section  is  devoted  to 
analysis  if  this  set  can  be  used  for  construction  of  conservation  laws  and  kinetic  equations 
similar  to  (57). 


5.  Model  of  multi-phase  flows  with  independent 

variables 


Now,  when  existence  of  a  single  thermodynamic  potential  is  obvious,  we  shall  try  to 
formulate  the  model,  which  would  involve  the  independent  variables  only.  We  select  the 
set  of  variables  (41)  at  z=1,...,m-1.  We  have  to  rewrite  the  specific  internal  energy 
applicable  to  the  present  case.  From  (47)  and  (63)  it  follows 


E  =  Z  CiEi=  Z  CiCi  +  CiV?-f2  =  Z'ciCi  +  Cnen+  (pV2)  ZciUi^  = 

=  Z 'Ci(ei  -  Cn)  +  Cn  +  (p2/2)ZCjl7i2  . 


(69) 


Turning  to  the  system  (57),  it  is  obvious,  that  the  mass  and  momentum  conservation 
equations  do  not  suffer  any  changes  associated  with  the  choice  of  independent  variables. 
The  energy  equation  contains  the  following  term  ZAjf^i,  which  changes  with  the  new 
choice  of  variables.  With  the  use  of  Z  Qi  =  0,  this  term  is  transformed  as  follows: 

ZAA  =  Z'AA  + An^2n  =  Z'AiOi-AnZ'^^i  =  Z'Qi(Ai-An)  .  (70) 

It  is  seen  that  use  of  the  function  Ai  -  An  is  preferable  in  the  present  case.  For  this  function 
to  be  involved,  we  have  to  modily  equations  for  the  specific  diffusion  velocity  by 
subtraction  of  the  last  one  from  the  rest  of  them  at  2=1,..., n-1.  We  introduce  a  new 
variable,  which  is  a  relative  diffusion  velocity  t  'i,  as  v'i=Vi-Vn-  Then  the  kinetic  equations 
for  follow  directly  from  (57): 

dpv\  ^  d(fyuv\  +  k\) . 

dt  dx 

Here  A'i  =  Ai  -  An  and  Wi  =  ^i  -  T^n.  Let  us  calculate  a  sum  of  the  diffusion  terms  in  (69), 
using  a  new  notation  for  the  relative  diffusion  velocity  and  (63): 

ZCiUi=  Z'CiUi+  CnVn  =  Z'CiUi+  (1  -  Z'Ci)  Un  =  Z'Ci(Ui-  Un)  +  t’n  =  Z'CiU'i+  Un  =  0  , 

Un  =  -Z'Cji;’j  , 

Ui=r’'i  +  yn  =y'i-Z'CjU’j  ,  (72) 
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Zcm^  =  Z'dVi^  +  CnVr?'  =  Z' Ci{v' i  -  Z’CjZ/j)2  +  (1  -  Z' Ci){Z' CjV' j)'^  = 

=  Z'Ciiv'iy  -  2  {Z'Ci  v'i)(Z'qv))  +  {Z' Ci)iZ' qv'jf  +  (2'qu’j)2  - 
-  (Z'ci)(2'qi;'j)2  =  Z' Ci{v' i)^  -  {Z' cp' j)^  . 

Thus,  the  set  of  independent  variables  (41)  is  transformed  into 

p,  S,  U,  Ci,  Bi,  v'i,  f=lr-  ■  ./tt-l.  (^3) 

Using  (71),  the  system  (57)  takes  the  following  form; 


dpc,  ^  d{puc.+Q^) 
dt  dx 


dpO,  ^  dpue, 
dt  dx 
dpv\  ^  d{puv;+A') 
dt  dx 


at  1=1,. .  .,n-l.  With  the  formula  (72)  the  specific  energy  (69)  is  written  out  as  foUows 

£  =  Z’Ci(ei-en)+e„+(p2/2)[Z'Ci(u’i)2-(Z’qi;'j)2]  .  (75) 

Finally,  we  have  to  adjust  the  thermodynamic  identity  (68)  to  the  variables  (73).  It  is  clear 
that  the  equations  T=Es  and  P  =  p'^Ep  preserve  their  form  because  the  transformation  does 
not  touch  scalar  variables.  Let  us  differentiate  the  specific  energy  (73)  over  the  variables  Ci, 
Bi,  and  u'i.  When  differentiating  over  Ciand  using  (63),  we  have 

£ci  ~  ~  Cn  Ci(£i)pi(pi)ci  "*■  Cn(Cn)pn(pn)cn  (Cn)ci  (P^/2)[(u  i)^  ~  2  (Z  Cp  j^V  J  — 


=  gi  -  Cn  +  pi{ei)pi  -  pn(en)pn  +  {p^/2)  v' i  [f'i  +  2Un]  = 

=  ei  -  gn  +  Pi/ Pi  -  Pn/pn  +  p^{vp  -  Vr?)/2  , 

which  is  exactly  Ai  -  An  =  A'i  from  (53).  Let  us  differentiate  £  over  Bi  and  use  (63): 
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Eqi  —  Ci(Ci)pi(pi)  0i  +  Cn(^n)pn(pn)  0n  (0n)0i  ~  P^/ P  P<^/ P  • 

Thus,  the  function  Hi  in  (59)  has  to  be  transfonned  into 

n'i  =  ni-nn.  (76) 

The  last  step  is  differentiation  of  E  in  (75)  over  v'i  with  the  use  of  (72): 

Evi  =  [  CiU'i  -  (Z'Cji’'j)-Ci]  =  p2ci  [u'i  -  Z'cjU’j]  =  fPCi  Vi  .  (77) 

It  is  seen  from  (43)  that  Ey'i  is  exactly  coincident  with  Qi. 

Thus,  the  thermodynamic  identity  can  be  easily  written  out  for  the  present  case: 

TdS  =  dE  +  Pdy -Z’A'idci-Z't2idu'i-Z’n’id0i  .  (78) 

Thus,  if  the  potential  E  is  given,  the  functions  involved  in  the  system  (74)  are  explicitly 
calculated  with  the  help  of  the  identity  (78). 

It  is  interesting  to  note  that  the  present  case  is  directly  reducible  to  the  case  of  the  two- 
phase  media  at  n=2  derived  in  the  third  section  because  w  =  v'i,  and  A  =  A'l  and  =  Q'l. 

The  identity  (78)  gives  us  the  rules  for  calculation  of  the  functions  entering  (74)  via  the 
potential  E(p,  S,  Ci,  9i,  v'i): 

T=Es,  P=p2Ep,  A'i  =  Eci  ,  Qi  =  Ev'i  ,  n'i=E0i.  (79) 

Expanding  the  energy  equation  in  (74)  and  applying  (79),  we  can  obtain  the  following 
entropy  evolution  equation: 

pT—  =  -z'A'm*  -  z'n;a),  -  z'q,'p;  ,  (80) 

dt 

here  d/df  =  d/dt  +  u  -d/dx  is  the  particle  derivative.  Using  (80),  the  correctness  requires  the 
following  entropy  nondecreasing  condition  to  be  satisfied: 

-Z'A'im‘i-Z'n'i<I)i-Z'QiW’i>0  .  (81) 

An  equivalent  condition  could  also  be  derived  from  the  overdetermined  system  (57),  (58): 


yOP— =  -ZA,.w* 
dt 


Zn,(D.  -103?, 


(82) 
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here  relations  Z  ^2;  =  0  and  Z  m'i  =  0  have  been  used.  It  is  interesting  to  note  that  for  the 
choice  of  kinetic  relations  suggested  in  [2]  as  wz‘i  =  0,  Oi  =  -  pOi/ x,  =  -  K^i ,  the  entropy 
nondecreasing  condition  is  satisfied  automatically.  Moreover,  the  Onsager  principle  [9]  is 
satisfied,  because  the  quadratic  form  pZ  (FIO^/x  +  kZ  (f2i)2,  appearing  for  this  case  as  the 
right-hand  side  in  (82),  is  symmetrical. 


6.  Hyperbolicity  of  the  two-phase  model 


In  this  section  we  analyse  eigenvalues  of  the  system  (36).  Hyperbolicity  of  the  prototype 
system  (27)  has  been  considered  in  [2]  and  the  eigenvalues  have  been  calculated  for  the 
case  of  absence  of  diffusivity  {lo  =  0).  To  analyse  the  system  (36)  we  shall  rewrite  it  in  the 
following  form,  replacing  the  energy  conservation  law  with  tiie  entropy  evolution 
equation  (similar  to  the  equation  (80)  for  the  case  of  two  phases): 


dp  du  - 

—  +  p —  =  0, 

dt  6x 


ox 

^P,Bp  ^  P,  8c  ^  P,  d0 
p  dx  p  dx  p  dx 
dc  dp  O  dc  n 
dt 

de 


+ 


P,  dv  Ps  dS 
’  -  +  -  ^  — 


p  dx  p  dx 
.0  dO  n„  dv 


:0, 


p  dx 
dt^^' 
dt 

dv  ^  dp 

rtf  rlv 


+  ■ 


p  dx  p  dx  p  dx 


=  0, 


dc  Kg  dO  A„  dv 


(83) 


here  d /dt  =  8 /dt  +  ii-d/dx  -  is  the  particle  derivative,  and  we  have  neglected  the  right-hand 
sides,  because  tliey  do  not  affect  ihe  characteristic  behaviour  of  the  system.  The  system  can 
in  general  be  written  out  in  the  following  matrix  form: 


dU 

dt 


+  A 


dU 

dx 


=  0, 


(84) 


where  U={p,u,  c,  6,  s,  v],  and  A  is  matrix  of  coefficients  of  the  system  (83).  As  a  result,  the 
characteristic  equation  of  the  system  det(A  -  \Z)  =  0  for  eigenvalues  X  is  specified  to  the 
following  characteristic  determinant 
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r  z 

p 

1  0 

1 

0 

0 

0  "1 

U- 

p,ip 

X 

1  pjp 

pjp 

pjp 

PJP 

0 

1  z' 

njp 

0 

^Jp 

0 

0 

i  0 

. 

z 

0 

0 

0 

0 

1  0 

1 

0 

z 

0 

0 

!  ^dp 

KIp 

KIp 

z'  J 

where  x  =  «  -  A  and  x'  =  X  =  Evc/p.  Expanding  (85),  the  characteristic  polynomial 

takes  the  following  form 


(  1  ^ 

A,Q, 

-P 

1  P  ) 

_ 

+^P,n,A,  - 

P  P  P 


XrP^K^X' — KPc^pX' 

P  P  P 


(86) 


=  0 


Combining  the  first  and  fifth  terms  in  the  square  parentheses  and  memorizing  the  zero 
double  root  x  =  0/  the  polynomial  for  the  rest  four  roots  is  reduced  to  the  following 


Let  us  introduce  the  following  denotations: 

A  =  Epc,  B  =  Epv,  ©i  =  Ecc/p,  ©2  =  Ew/p,  P  =  Pp.  (88) 

Then,  from  the  consequence  (34)  of  the  thermodynamic  identity  we  can  obtain 

Ap  =  A ,  Pc  =  p^A ,  S2p  =  B ,  Pv  =  p^B ,  Ac/p  =  ©1 ,  t2v/ p  =  ©2  .  (89) 

We  use  the  following  denotations  throughout  the  section  P  =  Pp  and  =  0i.©2.  It  is 
assumed  that  Pp,  ©1  and  ©2  are  positive.  That  gives  us  the  following  set  of  conditions  in 
terms  of  the  equation  of  state 


(p2Ep)p  >  0 ,  Ecc  >  0 ,  Ew  >  0  . 


(90) 


Then  the  characteristic  equation  takes  the  following  form: 
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F(x)-(X^ -1%'"  -0.  ■0j)+2pAB^’-^(A^0,  +B=0j=O.  (91) 

With  use  of  (9),  (24)  and  (88)  we  can  calculate  one  of  the  coefficients: 

A-B  =  pv{k^  -P  +  p^v^)  .  (92) 

It  can  be  straightforward  checked  out  that  the  polynomial  is  exactly  coincident  with  the 
one  derived  in  [2].  The  present  pol)momial  has  a  simpler  form  just  because  of  the  term 
reduction  associated  with  the  generalization  of  pressure  and  energy;  the  present 
generalization  contributes  in  additional  terms  related  to  the  density  derivatives  of  the 
diffusion  components  which  are  involved  in  the  characteristic  polynomial  of  [2]  as  a 
separate  term.  It  has  been  shown  in  [2]  that  roots  of  the  characteristic  polynomial  are  real 
in  the  vicinity  of  a  state  with  xu  =  0{v  =  0  for  the  present  model).  In  fact,  this  state  gives  B  = 
0  and  go  =  0,  resulting  in  the  following  bi-quadratic  equation 

(x2-Z2)(x2-k2)-pA202  =  O  . 

The  present  simpler  form  of  the  characteristic  polynomial  allows  us  to  assess  the 
hyperbolicity  of  the  model  in  a  wider  range  of  parameters.  To  study  the  polynomial  roots 
we  will  check  the  number  of  sign  changes.  Firstly,  we  rewrite  the  polynomial  in  a  form, 
which  is  more  convenient  for  the  analysis  employing  new  variable  z  =  X  2-  Expanding 

(91)  with  help  of  (92),  we  can  obtain  the  equation  for  z: 

P(z)  =  z4  -  2p-z2  +  2q-2  +  r  =  0  ,  (93) 


where 

p  =  (/c2+  /2+  eo2/2)/2  , 
q  =  (/c2 _ l2).eo/2  +  p^vik^ -P  +  p^v^)  , 

r  =  _  (fc2  +12.  eo2/4)-go2/4  +  {k^  -P  +  p'^v'^)-p^v-eo  -  p(A2©2  +  B2©i)  . 

It  is  obvious  that  the  polynomial  (93)  has  the  positive  sign  at  large  enough  negative  and 
positive  values  of  z.  If  we  show  that  the  polynomial  is  negative  at  two  certain  points  and 
has  the  positive  sign  at  a  point  between  the  two  then  it  means  that  the  polynomial  has  four 
real  roots.  We  choose  the  following  three  points: 

z_  =  -  pV2 ,  Z+  =  pi/2  ,  Zo  =  q/p  . 
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P(zo)  =  qVp^  +  r  . 


If  we  demand  that  P(z-)  <  0,  P(z+)  <  0  and  P(zo)  >  0  then  the  polynomial  (91)  will  have  four 
real  roots.  To  satisfy  these  conditions  it  is  sufficient  that 

r>0  ,  p2>r+2|q|  pV2  .  (94) 

When  the  denotations  in  (94)  having  been  decoded  with  (88)-(89),  these  conditions  (94) 
along  with  (90)  give  necessary  conditions  of  existence  of  4  roots  of  the  characteristic 
polynomial  (91).  When  the  generalized  specific  energy  E  is  specified,  the  conditions  (94) 
can  be  straightforward  checked  out.  The  values  z.,  z+,  and  zo  can  be  used  for  the  root 
separation.  The  estimate  (94)  is  not  the  best  possible  one;  for  instance,  choice  of  q/  (2p)  as  zo 
would  give  less  restrictive  assessment  for  die  coefficient  r  in  (94).  However,  the  conditions 
(94)  demonstrate  that  there  is  a  fairly  wide  range  of  parameters,  which  guarantee  existence 
of  real  roots  of  the  characteristic  polynomial. 


7.  Discussion  and  Conclusion 

A  model  for  multi-phase  flows  has  been  built  up  and  the  thermodynamic  identity  derived, 
allowing  us  to  close  the  model  with  a  single  thermodynamic  potential  -  a  generalized 
specific  energy  E.  In  doing  so,  a  link  between  the  generalized  pressure  and  energy, 
involving  the  diffusion  terms,  has  been  established  for  the  one-dimensional  case.  A 
condition  of  hyperbolicity  has  been  formulated  for  the  two-phase  model  for  a  range  of 
parameters  specified  by  the  equation  of  state. 

Probably  it  is  possible  to  generalize  a  three-dimensional  variant  of  the  model  [2]  to  the 
multi-phase  case  (in  fact,  after  the  manuscript  has  been  prepared,  the  author  received  an 
information  from  Prof.  E.  Romensky  that  a  generalization  of  the  prototype  model  [6]  has 
been  conducted  in  [10]);  however,  the  thermodynamic  identity  cannot  be  generalized  in 
the  same  manner  because  in  the  present  case  the  force-energy  link  is  scalar  and  this  fact 
has  been  essentially  employed  in  the  present  report;  whereas  in  the  three-dimensional  case 
this  link  is  of  essential  tensorial  nature. 

It  would  be  interesting  to  check  out  if  the  generalized  specific  energy  E  could  be  selected 
as  a  general  dependence  upon  the  specific  diffusion  velocity  Vi.  The  present  selection  as  the 
quadratic  dependence  reflects  consistency  with  the  inter-phase  behaviour  and  it  has  been 
used  for  the  design  of  the  model,  but  this  dependence  is  not  necessary  from  the  point  of 
view  of  diermodynamics  of  the  averaged  mixture.  In  that  case  this  dependence  could  be 
treated  as  a  specific  form  of  the  generalized  equation  of  state. 
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